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We study Bd and Bs mixing in unquenched lattice QCD employing the MILC collaboration 
gauge configurations that include u, d and s sea quarks based on the improved staggered quark 
(AsqTad) action and a highly improved gluon action. We implement the valence light quarks also 
with the AsqTad action and use the nonrelativistic NRQCD action for the valence b quark. We 
calculate hadronic matrix elements necessary for extracting CKM matrix elements from experimental 
measurements of mass differences AAId and AMs- We find ^ = /s, \/ Bb^ I fB^ \/ Bb^ = 1.258(33), 

/s^ \/ -Bbj = 216(15)MeV and /sj, \/ Bb^ ~ 266(18)MeV. We also update previous results for decay 
constants and obtain /s^ = 190(13)MeV, /s, = 231(15)MeV and /s./Zs^ = 1.226(26). The new 
lattice results lead to updated values for the ratio of CKM matrix elements \Vtd \ / \ Vts\ and for the 
Standard Model prediction for Br{Bs jj.'^fi^) with reduced errors. We determine |Vtd|/|Vfs| = 
0.214(1)(5) and Br{Bs /i+M") = 3.19(19) x 10"^ 



PACS numbers: 12.38.Gc, 13.20.Fc, 13.20.He 
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I. INTRODUCTION 

The mass differences AMg and AMd between the 
"heavy" and "light" mass eigenstates in the neutral B 
meson system have now been measured very accurately 
leading to the possibility of a precise determination of 
the ratio of two Cabibbo-Kobayashi-Maskawa (CKM) 
matrix elements |Vtd|/|Vts| P, 0, [1]. This ratio is an 
important ingredient in fixing one of the sides of the 
"Unitarity Triangle" , and hence plays a crucial role in 
consistency checks of the Standard Model. Reaching 
the goal of determining \Vtd \ / \ Vts\ from the experimen- 
tal AAfq's, however, requires theory input on hadronic 
matrix elements of certain four-fermion operators sand- 
wiched between the Bq and Bq states. Information on 
such hadronic matrix elements can only be obtained if 
one has control over the strong interactions, QCD, in the 
nonperturbative domain. 

In this article we use lattice QCD methods to calculate 
the hadronic matrix elements that appear in the Stan- 
dard Model to describe neutral B meson mixing. Simu- 
lations are carried out on unquenched configurations cre- 
ated by the MILC collaboration [4] . These configurations 
include effects from vacuum polarization due to three 
light quark flavors, up, down and strange {Nf = 2-1-1 
configurations, where Nf equals the number of sea quark 
flavors) . The up and down quark masses are set equal to 
each other. Table I lists the six different ensembles used, 
together with their characteristics such as the number 
of configurations, sea quark masses in lattice units, the 
valence quark masses employed for each ensemble, num- 
ber of time sources and the number of smearings for the 
b quark propagators. Information on the lattice spac- 



TABLE I: Details of configurations employed. Ntsrc is the 
number of time sources used per configuration and Nsm the 
number of smearings on the heavy propagator including the 
unsmeared local case. All quark masses are given in the MILC 
collaboration normalization convention with uo = (plaq.)^^'^. 
Errors in ri/a are estimated to be at the 0.5% level. 
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ing a is presented in terms of the ratio ri/a, where ri 
is obtained from the static potential and ri/a has been 
calculated by the MILC collaboration for each of their 
ensembles [51] . The bare b and s quark masses have been 
fixed already in previous simulations of the T and 
Kaon 7| systems. The MILC collaboration unquenched 
configurations have been created using the "fourth root" 
procedure to remove the four fold degeneracy of staggered 
fermions and some theoretical issues remain concerning 
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the validity of this procedure. Considerable progress has 
been made, however, in addressing this important issue 
and several recent reviews Q summarize our current 
understanding of the situation. In this work we assume 
that physical QCD is obtained in the continuum limit, as 
implied by existing evidence. 

In a previous article the HPQCD collaboration presented 
the first Nf = 2 + 1 unquenched results for Bg meson 
mixing parameters, based on simulations on two out of 
the above 6 MILC ensembles (sets C3 and C4) In 
the present work we broaden considerably the scope of 
our studies of B mixing phenomena. We generalize to 
include both Bd and Bs mixing and we use two sets of 
lattice spacings (the first 4 ensembles in Table I have 
a ~ 0.1 2fm and are called "coarse" whereas the last two 
have a 0.09fm and are refered to as "fine" lattices). 
Furthermore we now employ smeared operators for the 
Bq meson interpolating operators and even on those en- 
sembles used previously we have doubled the statistics, 
by going from two to four time sources. Unquenched 
lattice calculations by other groups exist in the litera- 
ture. Several years ago the JLQCD collaboration pub- 
lished Nf = 2 studies of Bj, and Bs mixing [ll| and 
the Fermilab/MILC collaboration has recently presented 
preliminary Nf = 2-1-1 results based on the same Asq- 
Tad light quarks as in the present article, however using 
a different action for the b quarks [l^, U^l • 

In the next section we summarize the formulas needed 
for analysis of B meson mixing phenomena. We intro- 
duce the relevant four-fermion operators and describe 
how their matrix elements are parameterized and how 
they can be related to the CKM matrix elements \ Vtd \ and 
|Vts|. We then discuss the lattice four-fermion operators 
used in the simulations and how they can be matched 
onto the operators in continuum QCD. In section III we 
describe our simulation data and the fitting procedures 
one must go through in order to extract the matrix ele- 
ments of interest. Section IV focuses on chiral and con- 
tinuum extrapolations and section V presents results for 



C = /b» Vb^s I /Sd a/SbT, /Sd V and /s^ Bb, to- 
gether with discussions of systematic errors. This section, 
section V, summarizes the main results of the present 
work for quantities most directly associated with B mix- 
ing analysis. As part of our simulations, however, we 
have also accumulated more data on B meson decay con- 
stants, /s^ and Jbs ■ Hence in section VI we update the 
results for these decay constants published previously in 
[l3 . [l5| . Section VII presents a summary of the current 
work and a discussion of future directions in our program. 

We conclude this introductory section with a comment 
on notation. The decay constants /s^ with q — d, s 
are defined in eq.(2) below and are used together with 
appropriate bag parameters to parameterize four-fermion 
operator matrix elements, /b^ can of course be identified 
with Jb, the decay constant of the charged B mesons, and 
/b can be measured through the latter meson's leptonic 



decays. The Bg meson, on the other hand, cannot decay 
leptonically via a single W boson and hence /s^ by itself 
is not a directly measurable quantity in the Standard 
Model. Although /b is the more physical quantity we 
use the notation /^^ throughout this article in order to 
facilitate uniform treatment of B^ and Bg mixing. 



II. RELEVANT FOUR FERMION OPERATORS 
AND MATCHING 

Neutral B meson mixing occurs at lowest order in the 
Standard Model through box diagrams involving the ex- 
change of two W bosons. These box diagrams can be 
well approximated by an effective Hamiltonian expressed 
in terms of four-fermion operators. More specifically, for 
calculations of AMq in QCD one is interested in the op- 
erator with [V-A] X [V-A] structure, 



0L = 



^i{v-A)m 



(1) 



where i and j are color indices and are summed over. 
The symbol q stands for cither the down or the strange 
quark. Working in the MS scheme, it is customary to 
parameterize the matrix element of OL between a Bq 



and a Bq state as. 



\MS' ^ ^ 



{OLy^ip) ^ {Bq\OL\Bqr^{y) ee -/^^ BbM^,- 

(2) 

Here Jb^ is the Bq meson decay constant and Bb^ its 
"bag parameter" . Factors such as | ensure that Bb, = 1 
in the "vacuum saturation" approximation. Given the 
definitions in ^ the Standard Model prediction for the 
mass difference is [16j], 



G\Ml 



W I 



VtqV:,\^i^iSo{xt)MBjlBB,. (3) 



67r2 



where Xt — /M^ depends on the top quark and the W 
boson masses mt and Mw, is a perturbative QCD cor- 
rection factor and So{xt) the Inami-Lim function. Bb^ is 
the renormalization group invariant bag parameter and 
at two- loop accuracy one has Bb /Bb = 1.539 for the 
present case. From eq.Q one sees that an experimen- 
tal measurement of AMq would yield directly the CKM 
matrix element combination iVt^Vj^p provided the quan- 
tity Jb^Bb^ is available. One also sees that the ratio 
|Vt(i|/|Vts| can be obtained from. 



\Vtd\ 



lAMdMB^ 
AM, Mb. 



^ fs, \/ Bb, 

/Sd \J BBa 



(4) 



The goal is to evaluate the hadronic matrix element in 
eq.([2]) using lattice QCD methods. Several steps are re- 
quired in going from what is actually simulated on the 
lattice to the MS scheme quantities appearing in the con- 
tinuum phenomenology formulas. One important step is 
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to relate four-fermion operators in continuum QCD to 
operators written in terms of the heavy and hght quark 
fields appearing in the lattice actions that we employ. 
Another crucial step will be to correct for the fact that 
simulations are carried out at nonzero lattice spacings 
and with light quark masses larger than the up or down 
quark masses in the real world. In the remainder of this 
section we address the first step, namely matching be- 
tween the continuum QCD operator OL and its counter- 
part in the effective lattice theory that we simulate. The 
other step of chiral and continuum extrapolations will be 
discussed in section IV. 

Our simulations are carried out using the improved stag- 
gered (AsqTad) quark action for the light quarks [13] and 
the nonrelativistic (NRQCD) action for the heavy quarks 
Matching through 0(q;^, Aqcd/M, as/(aM)) for 
the lattice action of this article was completed in refer- 
ence [l^, where M is the heavy quark mass. We refer the 
reader to that paper for details and just summarize the 
most important formulas here. In effective theories such 
as NRQCD one works separately with heavy quark fields 
that create heavy quarks (^'q) and with those that an- 
nihilate heavy antiquarks (^q)- The operator that con- 
tributes to Bq — Bq mixing at tree- level and that matches 
onto ID) at lowest order in 1/M has the form, 



%{V-A)^l ^'q{V-A)^1 (5) 



As is well known, even at lowest order in 1 /M there is a 
one-loop order mixing with another four-fermion opera- 
tor. 



^'QiS-P)^ 



(6) 



This is true both in NRQCD and in HQET. If one intro- 
duces an effective theory field. 



■eff 



Q 



*7 



(7) 



then and the QCD field are related by a Foldy- 
Wouthuysen-Tani (FWT) transformation. In particular, 



I 



—7 • V 

2M ' 



(8) 



where the V acts to the left. The FWT transformation 
determines the tree-level 1/M corrections to the four- 
fermion operators in the effective theory. For OL'^^^ they 
come in as. 



OL3I = — 



1 

2M 



¥q {V ~ A) -^q) (V% -^{V-A) ^q] 



Q 



(9) 



Taking these corrections into account one can work 
through 0{as,AQCD/M,as/{aM)) and finds the follow- 
ing matching relation, 

(0L)^(/.) = 

[1 + asPn] + as P12 {OS'^ff) + 



0{al,askQCD/M). 



(10) 



The matching coefficients pn, pi2, ^^^d C}"^ are listed 
(for /i = Mf,) in [l^. As explained there, the terms 
proportional to Q'^^ are needed to remove 0{as/{aM)) 
power law contributions in the matrix elements (OLjl). 



III. SIMULATION DATA AND FITTING 

The starting point for a lattice simulation determina- 
tion of (O), with 6 = OL^fi, OS'^ff or OLjl, is the 
calculation of the three-point correlator. 



.(4/) 

E 



C^p {tl,t2) - 

^ ,.,<i>^ (Xi,il) 0^(0) <^(f2,-i2)|0). 



(11) 



One works with dimensionless operators = a^O which 
are kept fixed at the origin of the lattice. <I>^ is an in- 
terpolating operator for the Bq meson of smearing type 
"a" , and spatial sums over xi and X2 ensure one is deal- 
ing with zero momentum Bq and Bq incoming and out- 
going states. The Bq meson is created at time — <2 and 
propagates to time slice where it mixes into a Bq me- 
son. The Bq meson then propagates further in time until 
it is annihilated at time ti. We have accumulated data 
for 1 < ti,t2 < Tjnax with Tmax — 24 On the coarse lat- 
tices and Tmax = 32 on the fine lattices. Given the well 
known properties of staggered light quarks, for fixed a, (3 
the three-point correlator must be fit to 

+ B"^ (-1)*' e"-^^'^*^"^^ g-Sfc(t2-i) 

j=0 k=Q 



j=Q k=0 



(12) 
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This ansatz allows for A'' non-oscillatory and N oscilla- 
tory contributions to the correlator (in practice we have 
worked with = N). Not all the amplitudes A^j^ etc. 



are independent due to symmetries, 
correlators are fit to, 



Similarly two-point 



N-l 



N-1 



(13) 



The relation between the amplitudes or the b" and 
the matrix elements of the previous section can be iden- 
tified as follows. 



A 



af3 
jk 



\^^\E,){E,\0^\E,){E,\<!>'il\Q) 
{2E,a3){2Eka^) ' 



(14) 



chi7dof 20 1.32 1.14 1.00 0.9 



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 
N 



The energy eigenstates in the numerator are taken to 
have conventional relativistic normalization and the fac- 
tors in the denominator are needed to make up the dif- 
ference between this continuum normalization and the 
one in the effective lattice theory. For the ground state 
contribution A^^ , and recalling that — a^O, one has, 



A 



00 _ 

00 ^ 



{B,\0\B,) 



|0) 



(2MbJ2 



(15) 



which includes the matrix element {Bq\0\Bq) that we are 
interested in. Similarly for the 2pt-functions one has, 



Using 



{B,Hl 



|0) 



(2MB,a3) 



(16) 



\njB,) 



one then has. 



Q/3 



{Bq\0\Bq) = 



2Mb, Aj^ 



(17) 



In order to assemble all the terms on the RHS of (iTOl) 
we have tried two approaches. In the first approach we 
did separate fits for each of the operators O = OL'^-^-^ , 
OS'^-^-f and OLjl and inserted their ground state ma- 
trix elements into (|10p . In the second approach we went 
through the analysis in the opposite order. Namely we 
first obtained the rcnormalized four-fermion operator at 
the three-point function level by forming the appropriate 
linear combinations of the C^'^-'^^'s, and then carried out 
fits to extract Aqq for the full rcnormalized three-point 
function. Consistent results were obtained from the two 
methods. For our final analysis we adopted the second 
approach which we found to be more convenient in prac- 
tice. 

Our smearings consist of Gaussian smearings of the heavy 
quark propagator at both source and sink. In addition to 



FIG. 1: Fit results for -W- versus the number of exponentials 

Nexp — N + N for one of the coarse ensembles, Set C2 with 
auoniyai = 0.04. 
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FIG. 2: Same as Fig.l for a fine ensemble, Set Fl with 
auonivai = 0.031. 



point sources and sinks we use Gaussians with widths, in 
units of the lattice spacing, of 2.0 and 6.0 for sets CI and 
C2 and one Gaussian each of width 6.0 for C3 and C4 and 
of width 8.0 for sets Fl and F2. To extract for our 
rcnormalized three-point function we carry out simulta- 
neous fits to an Nsrn X Nsm matrix of two-point correla- 
tors (eq. (fT3ll with a, /3 = 1, ...TV^™) (a = 1 corresponds to 
local, a = 2 to first Gaussian etc.) and to the rcnormal- 
ized three-point functions with a = f3. Bayesian fitting 
[20| methods are employed to enable these complicated 
fits with large numbers of exponentials, i.e. of fit parame- 

^1 — ^rnax 



ters. We fit to all data points within < t 



and i„ 



< f' for f 

— ''max '^^^ ''mm 



3, 



20 - 24 
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TABLE II: Fit results for fg^ ^ Mb^Bb, in units of and 

for the dimensionless ratio £< / a ,"' ~ V_^£f Errors 

in the last column are statistical + fitting errors. Those in 
the second and third columns include additional errors coming 
from the 0.5% uncertainty in ri/a. 



Set 




rV"fB,VMB,BB, 


V mZ 


CI 


1.430(21) 


1.193(27) 


1.199(29) 


C2 


1.442(16) 


1.248(35) 


1.155(33) 


C3 


1.382(21) 


1.179(21) 


1.172(23) 


C4 


1.413(18) 


1.263(22) 


1.119(22) 


Fl 
F2 


1.353(17) 
1.334(20) 


1.138(28) 
1.193(27) 


1.189(26) 
1.118(24) 



and Can = 13 15. We have used N^xp = N + N 
ranging between 4 to 16 and looked for consistency in fit 
results as the number of exponentials was increased. An 
example of fit results on one of the coarse lattices is shown 
in Fig.l. One sees that good and consistent results are 
obtained for 8 < N^xp < 16. When N^xp becomes very 
large (in the case of Fig.l > 16), errors tend to increase 
again indicating that our fit ansatz has become too com- 
plicated for the minimization routines to handle, given 
the amount of statistics that we have. Fig. 2 shows an 
example for one of the fine lattices. Here we find good 
results for 8 < N^xp < 14. In general we have relied on 
our N^xp = 8, 10 and 12 fits for all our ensembles. 

We summarize fit results in Table II. The dimensionful 

/ ~ —3/2 

quantities /b^ -w Mb^Bb^ are given in units of . Er- 
rors include both statistical plus fitting errors and errors 
coming from uncertaintiy in ri/a which we take to be 
~ 0.5%. Note that we have also gone to the renormal- 
ization group invariant bag parameter Bb^- 



IV. CHIRAL AND CONTINUUM 
EXTRAPOLATIONS 

The lattice data presented in Table II are for simula- 
tions with up and down quark masses m„ = vrid larger 
than in the real world and need to be extrapolated to the 
physical point. Reaching this physical point also involves 
taking the lattice spacing a ^ limit. We use stag- 
gered chiral perturbation theory (SChPT) [2l|, [H, ^ 
augmented by further general discretization correction 
terms to carry out the simultaneous chiral and contin- 
uum extrapolations. Continuum heavy meson chiral per- 
turbation for B and D mixing was developed in [23, HI] 
including for the partially quenched case. These formulas 
were generalized recently to next-to-leading order SChPT 
by Bernard, Laiho and Van de Water ^2^ and generously 
made available to us prior to publication. We use the fol- 



lowing fit ansatz, 

tTIb^^MbBb, = 

ci [1 + ^ A/g + C2 (2m/ + ms) ri + C3 niq ri] x 
[l + C4a,{a/n)^ +C5{a/n)% (18) 

A/g stands for the chiral log contributions and includes 
the staggered light quark action specific taste breaking 
terms. The factor of 1/2 comes about since Afq was cal- 
culated for the square, namely for fg Mb^Bb^- We use 
the notation m/ and for the sea up/ down and strange 
quark masses respectively, and (or to^s) for the va- 
lence quark masses. The second bracket parametrizes 
further discretization corrections that are expected to 
come in at O(asa^) and 0{a'^). We have also tried 
adding more analytic terms with higher powers of quark 
masses. 

A/g includes the coupling gB'Bir which has not been 
measured experimentally. However, based on Heavy 
Quark Effective Theory (HQET) arguments, gB'B-n is 
believed to be close to an analogous coupling gwD-n in 
the D meson system for which some experimental infor- 
mation is available. The latter coupling is estimated to 
be between 0.3 < go'Dir < 0.6 [23|. As we discuss be- 
low, we have carried out two types of fits, one where we 
did a whole sequence of fits with gB'Bir varying between 
< gB'BiT < 0.6 but where this coupling was kept fixed 
during each individual fit. In the second type of fit we let 
the coupling float and be one of the fit parameters. Both 
types of fits favored smaller values with gB'Bir ~ 0.1, 
however as long as gB'Bir < 0.5 fit results were quite 
insensitive to its exact value. 
For the ratio £,\/Mb, I^'^b, we use, 

[1 + i (Ao/g, - Ao/g) + ^ (Ai/g, - Ai/g) 

62 ijnqa - "ig) ri + 63 (rrigs - m^f- r^] x 

[1 + (64 a,(a/ri)2 + (a/ri)4)(mg, - m,) ri].(19) 

Here we have split up, 

A/g = Ao/g+5|.B.Ai/g (20) 

and then let g^. btx ~^ &i become one of the fit param- 
eters. In Fig. 3 we show a simultaneous fit to the six 
entries in the last column of Table II. The green and 
blue curves are the curves from this fit appropriate to 
the coarse and fine lattice data points respectively and 
the red curve is the "continuum" curve obtained by re- 
taining the fitted values for h\ , 62 and 63 and turning off 
the 64 and 65 correction terms plus the taste breaking 
contributions inside A/g and A/g^. One sees that within 
our statistical and fitting errors of ~ 2%, there is consis- 
tency between the three curves. In other words, we see 
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1.5 1 : I : I : I : I I I : I : I r 



♦ Coarse Lattice 
A Fine Lattice 
— Continuum 




1.5 1 : I : I : I : I I I : I : I r 



Coarse Lattice 
A Fine Lattice 
• Physical point 




oqI < I < \ < I < \ I \ < I < \ < I 

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 

r, m 

FIG. 3: Chiral and continuum extrapolation of the ratio 

^\/ '• Errors shown are statistical plus fitting errors. The 

physical point is at rimqs/27.4, where m^s is the valence 
strange quark mass. 



ngl < I < \ < I < \ I \ < I < \ < I 

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 

r, m 

1 q 

FIG. 4: Chiral and continuum extrapolation of 

3/2 

^1 fsdy Bb^Mb^- Errors on the data points are statisti- 
cal plus fitting errors combined with uncertainty in ri/a. The 
dashed curves correspond to the error bands about the central 
green and blue full lines. The physical point is at rimqs/27.4, 
where m^s is the valence strange quark mass. 



almost no statistically significant lattice spacing depen- 
dence in this ratio. At the physical point the difference 
between the green and blue curves is 1.8%, which reduces 
to 1.3% if the green curve is adjusted and corrected for 
having a sea strange quark mass on the coarse lattices 
that is about 20% too large. One might be surprised 
that the magenta curve lies below the blue curve. This 
comes about because the various discretization effects in- 
side {Afqs — Afq) and in the &4 & terms can have differ- 
ent signs and come in with different relative weights be- 
tween the coarse and fine lattices. All these effects come 
in at the ^ 0.5% or less level, and are hence too small to 
allow us to disentangle one from the other in a meaning- 
ful way. The fit shown in Fig. 3 has X^ug/dof = 0.54 [28j 
and gives gB*BTT = 0.14(47). 

Fig. 4 shows chiral & continuum extrapolation curves for 

^i^^/sd \/ Bb^Mb^ using the fit ansatz of eq. p5|) . Again 
the green and blue full curves are the fit curves for the 
coarse and fine lattice data respectively, and the dotted 
lines show the error bands around these central curves. 
Turning off the C4 and C5 contributions and the taste 
breaking terms inside A/^ leads to the red curve which 
can be followed down to the physical point. In contrast 

to the situation for the ratio ^, here, with Jb^ \J Bb^ ^^Bd i 
one finds a noticeable shift between the coarse and fine 
lattice points. The difference between the green and blue 
curves is a 5.5% effect. Going from the fine (blue) curve 
to the red continuum extrapolated curve is a 4% shift, 
which is also the size of the chiral & continuum extrap- 
olation error at the physical point. The fit in Fig. 4 has 
Xlug/dof = 0.99. 



3/2 / " 

Finally, in Fig. 5 we show results for r^' fB^y Bb^Mb,, 
where Xaug/d-of — 0.96 for the simultaneous fit to all the 
data points. Here the difference between the green and 
blue curves is a 6% effect and between the blue and red 
curve a 5.7% effect. These shifts are slightly larger than 
but similar to those for Bd in Fig. 4. In both cases the 

discretization effects we are seeing in r\^^ Jb^ iJ Bb^Mb^ 
are larger than the naive expectation of a leading cor- 
rection of ©(a^tts) which would be ~ 4% or ~ 2% on 
the coarse or fine lattices respectively. It was hence very 
important to have simulations results at more than one 
lattice spacing and carry out an explicit continuum ex- 
trapolation. Fortunately, for the important ratio ^ these 
discretization corrections cancel out to a large extent, as 
expected and as we have already verified in Fig. 3. 



V. MAIN RESULTS AND ERROR BUDGET 

Table III gives our error budget for the main uncertain- 
ties in the three quantities, /^^ Bbs : fsd \J BBd ^. 
We explain each entry in Table III in turn. 

• statistics and chiral extrapolations: These are the 
errors shown on the "physical points" in Figs. 3, 4 
and 5 and are outputs from our chiral & continuum 
extrapolation fits. 

• residual a} extrapolation error: It is necessary to 
list this error separately since the degree to which 
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FIG. 5: Same as Fig. 4 for r-j fs^y Bb^Mb^ versus 



nmf. 



TABLE III: Errors in % for /b, \/Bb, , fB^ V Bs^ and ^ 



source of error 




/sd V Bb^ 




stat. + chiral extrap. 


2.3 


4.1 


2.0 


residual extrap. 


3.0 


2.0 


0.3 


uncertainty 








r^^^ uncertainty 


2.3 


2.3 




Qb* Btt uncertainty 


1.0 


1.0 


1.0 


TTis and rrib tuning 


1.5 


1.0 


1.0 


operator matching 


4.0 


4.0 


0.7 


relativistic corr. 


2.5 


2.5 


0.4 


Total 


6.7 


7.1 


2.6 



the red curves in the above figures actually corre- 
spond to the true continuum limit depends on how 
well one has modelled discretization errors in our 
simulations. In other words one needs to assess the 
error in the fit ansatz for the continuum extrapola- 
tion (we assume the chiral extrapolation is han- 
dled suSiciently accurately by Staggered ChPT) 
and this turns out to be a nontrivial task. 
On the one hand the data appears to be consistent 
with the fit ansatze of eqs. (fT8|) and (fT9|) . We have 
tried adding further terms and found that fit results 
shifted by an amount less than, and in most cases 
much less than, the "statistical -I- chiral extrapola- 
tion" errors. Fig. 6 shows results for the chirally and 

3/2 _ 



continuum extrapolated r"]^ fs^^J Bb^Mb^ versus 
the number of parameters Npar in the fit ansatz. 
For Npar = 3 one has discretization corrections 
only through the Afg SxPT term and one sees 
that a good fit to the data points at the two lat- 
tice spacings cannot be obtained. However, once 
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Npar 



FIG. 6: 
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' Bb^Mb^ at the physical point versus 
Npar, the number of fit parameters. The numbers below 
the data points give Xaug/dof. For Npar = 6 we give results 
for two different types of term added, a term proportional to 
the quark mass squared (right data point) or a further dis- 
cretization correction cx {a/ri)^as (left data point). 



Npar > 3, good fits are achieved and results and 
their errors are stable with respect to changes in 
Npar. 

On the other hand we know that due to our use of 
the NRQCD action to describe the heavy b quark, 
coefficients such as C4 are in general complicated 
functions of (aM) (although C4 does include a con- 
stant piece coming from the light quarks). 
We have approached this complicated situation in 
the following way. We interpret the red "contin- 
uum" curves in Figs. 3, 4 and 5 as the curves one 
would get after taking care of all discretization er- 
rors coming from the light quark and the glue sec- 
tors. Then under "residual extrapolation un- 
certainty" one would include errors coming from 
the heavy quark action. The leading such error in 
our calculations is of 0{a'^as) multiplied by some 
function of (aM) which we initially take to be of 
0(1) leading to an additional uncertainty of ^ 2% 



in fsg y ■ This would correspond to a stan- 
dard power counting estimate of discretization er- 
rors where one takes coefficients of order one in 
higher order corrections. We have opted to be 
slightly more conservative in our power counting 
assessment and apply a factor of 1.5 rather than 
1.0 for those cases where the "physical" (the ma- 
genta) points deviate by more than one a from the 
fine lattice (blue) curve. By "ct" we mean here the 
"statistical -I- chiral extrapolation" errors. For ^ we 



have multiplied the error for /s^y Bb^ by a factor 
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of 



(m3-Tn„/d) 



1/6. 



Since the second row in Table III gives an assess- 
ment of the uncertainty in our a} extrapolation 
rather than an estimate of the full discretization 
error, we believe the procedure outlined here to fix 
it is a reasonably conservative one. 

• t^y"^ uncertainty: follows from the 1.5% error in 
current determinations of the physical value for r\ . 

• uncertainty in gB' b-k '■ we carried out fixed coupling 
chiral fits for the range 0.0 < gB'BTv < 0.6 and 
looked at the spread in the results at the physical 
point. For couplings larger than 0.6, /dof starts 
to deteriorate. 

• tuning of strange and bottom quark masses: The 
largest mistuning, which occurs in the sea strange 
quark mass mg on the coarse lattices, has been cor- 
rected for when calculating fit results at the phys- 
ical point and residual effects have been estimated 
by varying this adjusted value for TOs. Errors due 
to uncertainty in the valence strange quark mass 



have been assessed by comparing /s^ \J as one 
goes from valence quark mass mqs down to niq and 
errors coming from mistuning of m;, have been esti- 
mated from the 1/M dependence of decay constants 
studied in [l5|. 

• operator matching and relativistic corrections: 
These two sources of error are intimately inter- 
twined and again how to separate the two is not 
clear cut. As indicated in eq. pU]) . our matching for 
Jb Bb^ has been carried out up to correction of 
ofa^) and 0{asAQCD/M). In Table III we have 
listed the first correction under "operator match- 
ing" and the latter correction under "relativistic 
corrections" . And again the errors for ^ are re- 
duced by a factor of 1/6 relative to those for the 
two non ratio quantities. 

Using central values coming from the physical (red) 
points in the figures and the errors summarized in Ta- 
ble III, we can now present our main results. 



fBa\/BB^ 

and using ri = 0.321(5)/m 



1.258(25)(21), 



(21) 



Ib^^Bb^ = 266(6) (17) 



/b.V^b, =216(9)(12) 



0.321 

ri [Jm] 



0.321 

ri[/m] 
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MeV, (22) 



3/2 



MeV, 



(23) 



where the first error comes from statistics -I- chiral ex- 
trapolation and the second is the sum of all other sys- 
tematic errors added in quadrature. From the individual 



fB^y BBg, q=s or d, one obtains a ratio of 1.231 (58) (21) 
which is consistent with (|2ip however with larger errors. 



The result for fB^y Bb^ in ea. ip^ is consistent with 
but more accurate than our previously published value 
of 281(21)MeV [10]. 



VI. UPDATES ON fg^, fs, AND fg, / /s^ AND 
ESTIMATES OF BAG PARAMETERS 

The numerical simulations of two-point and three- 
point functions, such as in eqns. (fT3l) and (fTTj) . that en- 
abled us to extract the B-mixing parameters of the pre- 
vious section also provide information necessary to de- 
termine Bd and Bg meson decay constants /s^ and /b^ . 
Decay constants are defined through the matrix element 
of the heavy-light axial vector current between the Bq 
meson state and the hadronic vacuum. Using the tempo- 
ral component Aq and working in the heavy meson rest 
frame one has, 



{0\Ao\Bq)^MBjB, 



(24) 



Just as with the four-fermion operators of section II, 
matching is required between the heavy-light current in 
continuum QCD and currents made out of quark fields 
of the effective lattice theory. This matching has been 
carried out at the one- loop order for NRQCD / AsqTad 
currents in [2^ based on formalism developed in [30 |. 



;^o)^^ = (l + a.po)(Jr> + 

(l + a«Pi)(J^'^'^"') + agP2{4'^^'^') 
+ OialAlao/M'), (25) 

(i) 

The heavy-light currents Jq in the effective theory are 
defined as, 



7-(0)\ 



(0) 



J, 



y - 



r(2) _ 



2{aM) ^ ' 



(26) 
(27) 



^ *,7-V7oro*Q. (28) 



" 2 (aM) 

with Fq = 7570 and 

j{i),sub _ ^ j{Q) 

•Jq = "Jq ~ mo "^o 



(29) 



The matching coefficients pi and Qo are given in [29|. 
Note that the matching for the heavy-light current in- 
cludes contributions at 0{as ) and hence is more 



accurate than the matching in (jlOp for the four-fermion 
operator. 

We have evaluated the two-point functions. 



^(0|j(^)(xi,t)<^(f2,0)|0) (30) 
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for j = 0, 1, 2. We then calculate the renormalized cur- 
rent matrix element by forming the appropriate linear 
combination as dictated by the RHS of (|25p . This is done 
for both Bd and Bg. The next step is to fit the renor- 
malized two-point correlator using the ansatz of eq. (|13p . 
extract the relevant ground state matrix element and 
thereby obtain <i>^ = Jb^ ^Mb^- We do simultaneous 
fits to Bd and Bs correlators, so that ^d-, and the ra- 
tio <f>s/$(i are determined within the same fit. Fit results 

3/2 

are summarized in Table IV. For r-^ errors include the 
uncertainty in r-ija in addition to statistical and fitting 
errors. 

The rest of the analysis for $g and ^sl'^d is very similar 
to what was done for the four-fermion operator matrix 
elements in section IV. Chiral and continuum extrapola- 
tions are carried out using a fit ansatz of the form for 

r^x^^q and ([19]) for (f>s/$d. The only difference is that 
here A/, will involve the chiral logarithms appropriate 
for decay constants rather than for four-fermion oper- 
ators. Such contributions were calculated b y A ubin & 
Bernard using Staggered ChPT in reference [2^. Fig. 7, 
8 and 9 show chiral and continuum extrapolations for 
and with xlugldof = 1.00, 1.06 and 0.53 
respectively. 

Table V shows the error budget for Jb^ , /s^ and /s^ / /s^ , 
which is very similar to Table III for the mixing param- 
eters. The meaning of the different sources of error is 
as explained in section V. We have mentioned already 
that 0(as— Tj^) effects in the matching of the heavy- 
light current have been taken into account in our one- 
loop matching calculations ^29,]. Hence, these should not 
be included under "relativistic corrections" in Table V. 
However, there are still 0{as ^"^^ ) corrections to worry 
about in the NRQCD action. These would come from 
radiative corrections to the coefficient cb (often also de- 
noted C4) of the 2^ or • B term in the action. Although 
one-loop corrections to cb have not been calculated yet, 
one can nevertheless bound this coefficient nonperturba- 
tively by calculating the hyperfine, the B* — B, splitting 
and comparing with experiment. Preliminary results dis- 
cussed in fd\ indicate that cb is close to one and the en- 
tire effect would be at most a 10% correction to a —^7^ 
contribution. For the present calculations this means an 
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FIG. 7: Chiral and continuum extrapolation of the ratio 
c[>s/$d. The different curves and the physical point have same 
meanings as in Fig. 3. Here the red (continuum) curve is es- 
sentially on top of the blue (fine lattice) curve. 
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FIG. 8: Same as Fig.4 for r^^^^a = r^^^ fs^ ^/M^. 
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TABLE IV: Fit resuhs for <l>g = /s, \/Mb^ in units of 
and for the ratio ^a/^d- Errors are as described in Table II. 



Set 








CI 


1.261(12) 


1.085(14) 


1.162(14) 


C2 


1.246(11) 


1.073(14) 


1.162(12) 


C3 


1.236(12) 


1.071(14) 


1.155(14) 


C4 


1.248(16) 


1.128(17) 


1.107(20) 


Fl 


1.175(13) 


0.990(22) 


1.188(20) 


F2 


1.180(13) 


1.047(16) 


1.120(11) 



source of error 






fBjfB, 


stat. + chiral extrap. 


2.2 


3.5 


1.6 


residual extrap. 


3.0 


3.0 


0.5 


uncertainty 
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uncertamty 


2.3 


2.3 




Qb'Btv uncertainty 


1.0 


1.0 


0.3 


TUs and nil, tuning 


1.5 


1.0 


1.0 


operator matching 


4.0 


4.0 


0.7 


relativistic corr. 


1.0 


1.0 


0.2 


Total 


6.3 


6.7 


2.1 
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FIG. 9: Same as Fig.5 for ri^'^<!>s = rl^^fs, ^JMb,. 



uncertainty of order 1% in fs^ and a much smaller one 
for 1^. 

The final numbers for the decay constants including all 
errors added in quadrature become, 
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190(13) 



/b^, = 231(15) 



1.226(26), 



0.321 

ri [/to] 



0.321 
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MeV, 



3/2 



MeV. 
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FIG. 10; Chiral and continuum extrapolation of the bag 
parameter Bb,- In order to avoid clutter, we show an error 
band only for the fit to the coarse lattice data. 



and V are either irrelevant or cancel to a large extent 
in -B_B^. For instance, one sees from Fig. 10 that there is 
little evidence for discretization errors in the bag parame- 
ter, although /b^ ^^Mb^Bb, and /b, \/Mb, individually 
do have noticeable lattice spacing dependence. Our final 
result for the Bs meson bag parameter is Bb^ = 1-33(6), 
or BB,{Mb) — 0.86(4). One sees that direct extraction 
of the bag parameter followed by a continuum/chiral ex- 
trapolation can reduce errors significantly over our first 
approach of extrapolating Jb, \/Mb,Bb, and /s^ y/ M b^ 
first and then taking the ratio. Although we were success- 
ful in applying the second method for Bb^ , unfortunately 
it has not been possible to get stable extractions of Bb^ 
for all of our six ensembles. Our best estimate for the 
Bd meson bag parameter is obtained by taking the two 
ratios, fB.y^B^/fB^y^B^and /bJIb^, from eqs.(21) 
and (31) and combining their ratio with our most accu- 
rate BBAMb). This leads to BbJBb^ = 1.05(7) and 
Bb, (Mb) = 0.82(7). 



These results for Jb^ are consistent with but about 
one a lower than the values /r . = 216(22)MeV and 
Jb^ = 260(29)MeV given in [lll3- The main differ- 
ence between the analysis carried out here and in (l3 ] 
is that in the latter case chiral extrapolations were done 
based only on coarse lattice data and furthermore no at- 
tempt was made to extrapolate explicitly to the contin- 
uum limit. The new result for the ratio in (|3ip is similarly 
consistent with our previous /_b^//b^ — 1.20(3)(1) p^ . 

The Bq mixing simulations can also be used to deter- 
mine the bag parameters Bb^- From the separate fi- Our main results, namely values for ^ = 
nal results for Jb^ y^Bq S''^^ ^^r Jb^ one finds Bb, (/i = 
A/b) = 0.86(6) and BbAM^) = 0.84(10). We have also 
attempted to extract the bag parameters directly from 
simultaneous fits to 3-point and 2-point correlators for 
each ensemble separately. The results are shown for Bb^ 
in Fig. 10 and extrapolated to the physical point, where 
one finds Bb, — 1.33(5). We add to this 3.8% statistical 
-|- chiral extrapolation error additional 2.5% systematic 
errors. Many of the systematic errors listed in Tables III 



VII. SUMMARY 

We have completed the first Nf = 2 + 1 unqucnchcd 
study of Bs and Bd mixing phenomena in Lattice QCD. 



Jb, y Bb, and /^^ y Bb^ , are given in eqns.lUT]), and 
respectively. Combining the lattice result for ^ with 
the experimentally measured mass differences AAf^ = 
0.507±0.005ps-i d and AMs = 17.77±0.10±0.07ps-i 
[3] leads to. 



\Vtd\ 

\Vts\ 



0.214(1)(5) 



(34) 
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where the first error is experimental and the second 
from the lattice calculation presented here. This is the 
first time this ratio of CKM matrix elements has been 
determined while incorporating a fully self consistent 
Nf = 2 + 1 calculation of ^. 

In addition to giving mixing parameter results, in this ar- 
ticle we have also updated values for decay constants fs^ 
and /b, and their ratio in section VI. fs, appears in the 
Standard Model prediction for the branching fraction for 
Bs dHj ^ process sensitive to new physics. The 

most accurate result for this branching fraction comes 
from taking a ratio with AMg [13], which gives a result 
inversel y p roportional to Bb, ■ Updating the parameters 
used in [33| for t{Bs) and mt [H, and including our result 
for Bb, of 1.33(6), gives 

Br{Bs = 3.19(19) X 10"^ (35) 

improving on the previous accuracy available. The 
largest contribution to the error on the branching fraction 
comes from the error on Bb^ followed by the uncertainty 
in t{B,). 

The calculations presented here can be improved in sev- 
eral ways. Foremost among the improvements planned 
for the future is to carry out simulations at other, finer, 
lattice spacings. Having results at more than two lattice 
spacings will help considerably in reducing the "statis- 
tical -I- chiral extrapolation" and the "residual ex- 
trapolation" uncertainties in Tables III & V. They would 
also contribute to constraining the value of qb' b-k so that 
this source of error can then be ignored. Hence, one can 
expect lattice results for ^ (and also for JbJ fBa) with 
accuracy of ~ 1% in the not too distant future. Im- 
provement for dimensionful quantities such as Jb^ \J Bq 

3/2 

will also require reducing the "rj^' " and the "operator 
matching" errors. The HPQCD collaboration is currently 
engaged in projects aimed at fixing the physical value of 
Ti [33] with higher precision than in the past. We are 
also exploring nonperturbative methods for carrying out 
operator matching in heavy-light systems. At least for 
heavy- light currents, methods recently applied to accu- 
rate determinations of heavy quark masses, which involve 
moments of current correlators and very high order con- 
tinuum QCD perturbation theory (3^ 1 , look promising for 



nonperturbative determinations of Z-factors. More work 
will be necessary to see whether such methods can be 
generalized to four-fermion operators. It is possible one 
can take advantage of the fact that a major contribution 
to matching of four-fermion operators comes from dia- 
grams involving radiative corrections to just one of the 
bilinears within the four-fermion operator, in other words 
corrections that are identical to a heavy-light current ra- 
diative correction. This has been noted already in the 
one-loop calculations of [l^ . With several of these im- 
provements in place better than ~ 5% accuracy should 

be achievable for fs^ \J Bb^ ■ 

Another worthwhile direction for future investigations 
would be to calculate hadronic matrix elements of fur- 
ther AB = 2 four-fermion operators, beyond the two, OL 
and OS, studied here. As is well known, there are five 
suc h op erators usually denoted Qi, with i = 1,2,3,4,5 
[ll,[3l[3i|. In this notation OL = Q1 and OS = Q2. In 
this article we have focused on Ql and Q2 since only they 
are relevant for the mass difference AAIq in the Standard 
Model. The operator Q3 would come in for calculations 
of the width difference AF^ [s^l- Although we have al- 
ready accumulated simulation data for {Q3) we will post- 
pone analysis for a future publication where we also plan 
to have results for ((54) and ((35). In [19| the necessary 
matching at one loop has already been completed for all 
five four-fermion operators. The two hadronic matrix 
elements ((^4) and {Q5) do not appear in the Standard 
Model but are of interest in several Supersymmetric Mod- 
els. To date only quenched lattice results exist for all 
five four-fermion operators (36j . It will be important for 
Beyond the Standard Model studies to generalize those 
results to unquenched calculations. 
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